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Abstract 

In a previous investigation, a model of three-body motion was devel- 
oped which included the effects of gravitational radiation reaction. The 
aim was to describe the motion of a relativistic binary pulsar that is per- 
turbed by a third mass and look for resonances between the binary and 
third mass orbits. Numerical integration of an equation of relative motion 
that approximates the binary gives evidence of such resonances. These 
(m : n) resonances are defined for the present purposes by the resonance 
condition, muj = 2nQ, where m and n are relatively prime integers and 
u> and fi are the angular frequencies of the binary orbit and third mass 
orbit (around the center-of-mass of the binary), respectively. The reso- 
nance condition consequently fixes a value for the semimajor axis a of 
the binary orbit for the duration of the resonance because of the Kepler 
relationship to = a~ 3//2 . This paper outlines a method of averaging de- 
veloped by Chicone, Mashhoon, and Retzloff which renders a nonlinear 
system that undergoes resonance capture into a mathematically amenable 
form. This method is applied to the present system and one arrives at an 
analytical solution that describes the average motion during resonance. 
Furthermore, prominent features of the full nonlinear system, such as the 
frequency of oscillation and antidamping, accord with their analytically 
derived formulae. 
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1 INTRODUCTION 

A model was established in a previous paper, 'Gravitational Radiation Reac- 
tion and the Three Body Problem' (Wardell 2002), in which the motion of a 
binary system was studied. This binary system was subject to gravitational 
radiation damping and the gravitational influence of a third mass. The equa- 
tions of motion for the relative motion of the binary were derived with certain 
approximations that would highlight the effects under investigation and also 
make analysis easier. For example, the motion of the three masses is taken to 
be planar and the center-of-mass of the binary moves in a fixed circular orbit 
around the third mass. Furthermore, the distance of the third mass from the 
binary's center-of-mass is taken to be substantially larger than the size of the 
relative orbit of the binary. The masses are considered to be point masses. After 
a scaling transformation that makes the variables dimensionless, one arrives at 
the following form of the equation of motion in Cartesian coordinates: 

— = ---eK^-8W. (1) 

The variable r l represents the relative orbit, Kij(t) is the tidal matrix that re- 
tains the information about the tidal interaction between the third mass and the 
binary system, and R l is the radiation reaction term that expresses the radiation 
reaction force to desired order after iterative reduction. To avoid 'runaway so- 
lutions' that can arise as the result of the radiation reaction perturbation which 
involves a fifth time derivative, one can apply the method of iterative reduction. 
This gives rise to a set of equations that is second order and resembles an equa- 
tion of motion in Newtonian mechanics (Chicone et al. 2001). One recovers the 
differential equation for the appropriate Kepler problem if e = 6 = 0. 

Numerical analysis of the system has revealed that given appropriate initial 
conditions, one arrives at a result which shows resonance behavior in the relative 
orbit. When the graph of the Delaunay variable L, where L = a 1 / 2 such that a is 
the semimajor axis of the osculating ellipse of the relative orbit, is plotted versus 
time one sees that the trend of semimajor axis decay can temporarily stop on 
average at a resonance. An oscillation occurs around a fixed average value for 
L for the duration of this resonance. This resonance capture is indicative of an 
average balance of energy that leaves the binary system by way of gravitational 
waves and enters the system because of the tidal gravitational influence of the 
third mass. 

This paper concerns itself with the behavior of the system when a resonance 
occurs; that is, when the resonance condition muj = nil' is satisfied, where m 
and n are relatively prime integers and u> and f2' are the angular frequencies 
of the relative motion and tidal perturbation respectively. It turns out that 
the tidal perturbation frequency relates to the fixed third-body frequency as 
O' = 20, where ft is the frequency of the fixed third-body motion. 

An averaging method that was developed for the purposes of studying the 
effect of external incident gravitational waves on a binary system, can be applied 
to the equations of motion to analyze the behavior near a resonance (Chicone, 
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Mashhoon & Retzloff 1997). The resultant averaged system ol equations re- 
tains the 'slow' variables which are left over after the system is averaged over a 
'fast' variable. The averaged set of nonlinear equations gives rise to a solvable 
system, whose solution approximates the actual solution for sufficiently small 
perturbation parameters over a certain time scale. 

2 NUMERICAL RESULTS 

One arrives at the following system of equations upon converting the second- 
order Cartesian equations into a system of first-order equations in polar coor- 
dinates: 



r = P r , (2) 




P 2 1 1 P V) P 2 

P r = ^-- 2 + - 2 er[l + 3 cos (n't 29)] + 6^ + 24P 2 + 144-J-), (4) 

3 Pa 94 P 2 

P e = ?-er 2 S in(n't-2d) + 5^(— + l08P? -72-^). (5) 

The parameters e and 5 denote the tidal and gravitational damping perturbation 
strengths, respectively. One can in principle generate a solution to this system of 
differential equations by specifying a set of initial conditions that represent the 
physical state of the system. The system of equations also includes parameters 
e, <5, and fl that characterize the strength of the perturbations and the angular 
frequency of the tidal perturbation. These parameters can be picked to reflect 
a specific physical model. In fact, a simplification has been made which, in 
effect, makes e an independent small parameter. The original derivation reveals 
that the tidal perturbation amplitude is directly proportional to fi' 2 , which 
may not be small enough in many cases. A generalization is made in this 
model to accommodate the requirements of the averaging theorem, particularly 
that the perturbation amplitude be sufficiently small. A series of numerical 
experiments can be performed in which the system evolves under the influence 
of the prescribed perturbations and initial state. This process aims to find 
occurrences of resonance capture in phase space and parameter space of the 
dynamical system. 

To accord with the subsequent mathematical analysis, the Delaunay variable 
L is graphed versus time in the numerical work. The variable L is related to 
the semimajor axis a of the Keplerian ellipse as L — a 1 / 2 . Therefore it is 
inversely proportional to the energy of the elliptical orbit as L = (—2E)~ 1 / 2 . 
The standard signature of the resonance capture in terms of L versus t is an 
oscillating interval about an average value L„, with an envelope that generally 
increases with time before 'falling out' of the resonance. The average value is 
related to the order of the resonance. A sustained average orbital energy follows 
directly from the average L. 
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Numerical experiments are conducted with the aim of finding resonances of 
different orders. In this investigation, numerical evidence of (1 : 1), (2 : 1), and 
(3:1) resonances is found. See figures 1-3 and their captions which show the 
parameters and initial conditions of the systems. It is not difficult to search 
for resonances in the parameter space of the system; in particular, the (1:1) 
resonance described here is different from the one represented in a previous 
paper (Wardell 2002). A cursory look at the graphs will show that there is 
an increasing degree of structure, especially with regards to the (3 : 1) graph. 
One can see that the (3 : 1) resonance contains a 'dense' orbit with chaotic 
characteristics. That is, the features of the graph such as the oscillations are 
apparently irregular as is the amplitude. In contrast, the (1 : 1) and (2 : 1) 
orbits are much more regular in their oscillations. 

The graphs exhibit similar qualitative behavior to which a measurement can 
be applied. The oscillation during the resonance capture gives rise to an average 
frequency. Let this be to m . Moreover, the amplitude of the envelope of the 
oscillation changes with time and yields another measurable quantity. Namely, 
the ratio between the amplitudes of the envelope at two different times is a 
simple measure of the change of the envelope that surrounds the oscillatory orbit. 
Let this be R m . The purpose of this paper is to provide analytic expressions 
for u) m and R m and an approximate theoretical description of the binary orbit 
while in resonance. 

3 EQUATIONS OF MOTION IN DELAUNAY 
VARIABLES 

If one neglects the damping terms in the equations of motion (2)-(5), one arrives 
at a Hamiltonian system that can be transformed canonically into other coordi- 
nate systems. The coordinate transformation from polar to Delaunay variables 
is beneficial because Delaunay variables are action-angle variables adapted to 
the osculating ellipse. The Hamiltonian of the unperturbed system depends only 
on the actions or the canonical momenta. This results in the momenta being 
constants of the motion and the generalized coordinates being angle variables. 
The inclusion of perturbations causes deviations from these constants — the 
otherwise constant variables evolve with time. The close relationship that the 
Delaunay variables have with the elliptical elements conveniently connect the 
trajectory on the manifold with respect to the Delaunay variables with the os- 
culating ellipse description of motion that is associated with elliptical elements. 
In other words, one has at each point of time a set of coordinates that will 
follow the Keplerian trajectory if the external perturbations were removed. The 
orbit of the system is aptly described by an evolving ellipse, whose descriptive 
coordinates change in time. For the planar system under consideration here, 
the osculating ellipse can be characterized by its semimajor axis a, eccentricity 
e, eccentric anomaly u, and true anomaly v (Danby 1988). The Delaunay ele- 
ments are then defined by L = a 1 / 2 , G = Pg = L(l — e 2 ) 1 / 2 , t = u — esinu, and 
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Figure 1: This is a (1:1) resonance. The parameters are e = 5 x 10 5 , 
S/e = 10~ 3 , and 0' = 0.46. The initial conditions at t = are (r, 6, P r , P e ) = 
(1,0.5,0.63569,1). 
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Figure 2: This is a (2:1) resonance. The parameters are e = 5 x 10 6 , 
S/e = 1CT 3 , and O' = 0.30. The initial conditions at t = are (r,6,P r ,P e ) = 
(1,0.5,0.847165,1). 
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Figure 3: This is a (3:1) resonance. The parameters are e = 10 6 , 6/e — 
1(T 3 , and f2' = 0.10. The initial conditions at t = arc (r, 0, P r , Pq) = 
(1,0.5,0.72059,1). 
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g = 9 — v. Here L and G are action variables and the mean anomaly I and g 
are the corresponding angle variables. 

The equations of motion in Delaunay variables take the following general 
form 

L=-e^+eAf L , (6) 
G = -e^ + eAf G , (7) 

9 = e^f+eA fg , (9) 

where A = 5/e. This definition is made so that one can fix A and have only 
one small parameter e that is needed for the remaining analysis. The averaging 
theorem relies on the system having one perturbation parameter. 

The Hamiltonian associated with the binary system plus the external per- 
turbation due to the third mass can be expressed as 

H = Hq + eHext, (10) 

where the unperturbed part is 

"o = -^, (ID 

and the tidal perturbation is 

U ext = A{L,GJ) + B{L,G,e,g)cosQt + C(L,GJ,g)smQt. (12) 

The functions A, B, and C are given by 

1., 3. G 2 . cosu£ T , .. , . 

A=- 1 L\l + -(l-—)+J2^^Mve)), (13) 

V— 1 

15 G 2 3 °° 

B = -— i 4 (l - - 7 )cos2 5 - - L A [A v cos 2g cos vl — B u sin 2g sin vl) , (14) 

U—l 

15 G 2 3 °° 

C = — -L 4 (l - — ) sin 2^ - -L 4 sin 2# cos ^ + ^ cos 2# sin vl)), (15) 

o -L 4 — ( 

v—1 

where 

A, = 7^[2^e(l - e 2 ) J>e) - (2 - e 2 )J v {ve)], (16) 
B„ = — ^(1 - e 2 )i [eJ>e) v(l e 2 )J v {ve)]. (17) 
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The symbol J„ represents a Bessel function of the first kind of order v and 
e = (1 - G 2 /i 2 ) 1 / 2 is the eccentricity. 

The transformation for the general case which includes damping uses the 
following relationship which accounts for the non-conservative nature of the 
damping force: 

1 dE „ 

-— = F v, 18 

(i at 

where /x is the reduced mass. One can use this relation to find how the damping 
terms transform with the new coordinates. This result was derived in a previous 
study where the damping terms are the same as in this case (Chicone et al. 
1997). The terms are: 



Lr 3 



[1- 



16 L 2 



,20 



+ (^-L 2 - ^-G 2 ) 



17^L 2 50 L 4 G 2 25L 4 G 4 n 



+ 



fa = - 



18G 



[1- 



20 L 2 5L 2 G 2 . 



9 r 



+ 



(19) 
(20) 



2sinw r „ 2 l /wo -j ^^1 „J0 r , 27^ 2N G 2 25 L 2 G 4 ^L 2 G 6 , 

ft = ~ r37r"2 [4e ~r--(73G — 40.L ) — 2(—L ——G )—^ 3 — 1-25 — —], 

eL^Gr 1 3 r 3 2 r l 3 r A r 4 

(21) 

2sinw r , 2 80 2x l 25i 2 G 2 i 2 G 4 n , N 

h - -^3 HI + - y L 2 )- - y - 3 - + 25-^]. (22) 

So one has here all of the ingredients to produce the equations of motion in 
Delaunay variables. The system is formally suited for the averaging process. 



4 PARTIAL AVERAGING NEAR RESONANCE 
4.1 FORMALISM 

When the perturbation is included, the action variables L and G as well as the 
angle variable g change in time along with the angle I. Their time derivatives 
in the equations of motion (6)-(9) are equal to an expression of the order of the 
small parameter e; therefore, they evolve slowly over time. The angle variable I 
maintains its leading term which is of order unity — so this describes fast motion 
in relation to the other variables. In essence, one is interested in the evolution 
of the slow variables L, G, and g averaged over the fast motion of £. The motion 
of I for the unperturbed system is just the 'mean angle' for the relative orbit. 
The variable L is related to the energy of the orbit by E — —1/(2L 2 ), G is the 
angular momentum, and g is the argument of the periastron of the osculating 
ellipse. Averaging over the fast angle brings to light an averaged solution of the 
slowly evolving variables, particularly observable variables such as L. 
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To advance an analytical description of the modelled system, one can take 
advantage of the qualitative features of the ODE — namely, the fact that there 
is one 'fast' angle which evolves quickly with respect to the 'slow' variables. 
The 'slow' angle variable can be considered to be effectively an action variable. 
This is key because to meet the conditions of the averaging theorem outlined 
below, one needs a system with only one fast angle variable. If one makes use 
of the fact that there is only one 'fast moving' angle in the system of interest, 
the averaging theorem can be applied. Employing the averaging theorem, one 
arrives at an averaged set of differential equations that are solvable. Over the 
correct timescale, one has an approximate analytical solution to the motion of 
the binary system of interest. 

To illustrate how averaging is applied using the averaging theorem consider 
the following representation of a system of equations cast into action-angle form: 

i = ef(I,6), (23) 
9 = n(I) + eG(I,d), (24) 

where 



/,fer, (25) 
e,n,ge$t. (26) 

The symbol Sft represents the set of real numbers and 5ft" represents an n- 
dimcnsional real space. The terms with the coefficient e are perturbations. 
Furthermore, the evolution of the action variables I is slow in time with respect 
to the angle variable 9 that is fast moving. So one has a case where averaging 
can be applied. One can average over the angle 9 and create a new differential 
equation so that to first-order one has a new system in terms of J: 

j = e±- JJ F{J,6)de. (27) 

The new averaged variable J coincides with the actual variable I at t = 0: 

J(0) = 1(0) = I . (28) 

The Averaging Theorem states that one can create an averaged system whose 
solution will closely follow the actual system's solution. That is the trajectory of 
the new averaged variable J(t) will stay close to the true trajectory I(t) within a 
distance of order e. This will be valid within a timescale of e _1 . So the following 
inequalities describe the averaged system where C and C are constants: 

\I(t)-J(t)\<Ce, (29) 
C 

< t < — . (30) 

Inspection of the differential equation in the Delaunay variables reveals that 
the 'fast' variable is £ while the others are 'slow'. It is of interest to study the 



10 



behavior of the system's motion near a resonance. At a resonance, as can be 
seen in the numerical case, the value for L stays fixed on average. For a fixed L* 
associated with a particular resonance, one arrives at the following resonance 
condition in terms of L„: 

= nQ! (31) 
where lo = L~ 3 . The unperturbed equation for I yields the solution 

£(t) = ^t + £ , (32) 

where l Q is an integration constant. 

To highlight the motion near resonance, one can consider the deviations 
that occur near the resonance. The averaging that one performs at resonance, 
is called partial averaging. To investigate the behavior of the binary near the 
resonance manifold G,£,g) one makes the transformation 

L = L„ + e 1/2 L>, £ = j^ + <P- ( 33 ) 

This introduces the variable D which represents the 'deviation' from the reso- 
nance value of L* and the variable ip which is related to a deviation from the 
unperturbed mean angle I. The transformation involves e 1 / 2 before the D vari- 
able and unity before the ip variable so that the transformed ODE will contain 
the differential equations in D and if with leading terms to the same order. The 
leading order becomes the small parameter for the averaged system, that is e 1 / 2 . 
This gives the proper form for averaging. 

The derivation involves an expansion around L*, therefore it is possible to 
mathematically generate any number of terms in the expansion. Since the model 
equations that reflect the physics under investigation are given to highest order 
e, one does not benefit from a physics point of view to expand beyond e in the 
transformation. 

After the transformation is made one obtains the following expression for 
the transformed system at resonance (Chicone et al., 1997): 

D = -e 1/2 F u - eDF 12 + 0(e 3/2 ), (34) 
G = -eF 22 + 0(e 3 / 2 ), (35) 

V = -t 1/2 ^ D + <^, d2 + ^32) + O(e^), (36) 

g = eF 42 + 0(e 3 / 2 ). (37) 

Here Fij arc defined by: 



Fn := ^(L„G, —t + ip, g, t) - Af L (L*, G, —t + ip, g), (38) 



li 



F 12 :=^(L*,G,—t + ^g,t), (39) 



F 22 := ^N(L., G, — t + <p, g, t) A/ G (L„, G, —t + <p,g), (40) 
og m m 

F 32 := ^(L., G, —t + <p, g, t) + A/,(L*, G, —t + <p,g), (41) 
oL m m 

F 42 := (L., G, — i + V) .g, t) + A/ g (L*,G, — t + 5 ). (42) 

The next step is to derive the actual second order partially averaged equa- 
tions for the system. To do this, one can make use of a special coordinate change 
known as an averaging transformation. The purpose of such a transformation 
is to render the system by way of a coordinate change into a form that is aver- 
aged to first order. Since averaging to second order is desired here, one needs 
to average the entire transformed system and drop terms of order higher than 
second. 

To apply the averaging transformation one must make the following defini- 
tions: 



- ft' f 2 ^ _ nft' 

Fij := / Fij(G, s + ip,g lS )ds, 



(43) 



such that 



nil' — 
\(G,<p,g,t) ~Fn(G,—t + <p,g,t)-F u , (44) 

A(G,(p,g,t) := J \(G,<p, g, s)ds, (45) 

-Sap 

/ " A(G,<p,g,s)ds = 0. (46) 
Jo 

That is, the average of A vanishes. The averaging transformation becomes: 

D = D-e 1 / 2 A(G,<p,g,t), G = G, <p = <p, g = g. (47) 

This transformation has the property that its average becomes the identity 
transformation. The averaging transformation changes equations (23)-(26) into 
the following: 

D = -e^Fn eD(F 12 + + O(e^), (48) 

G = -eF 22 + 0(e 3 / 2 ), (49) 
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* = -' 1/2 h b + t{ ij b2 + F " 2 + h K) + 0(e3/2) ' (50) 

^eF 42 +0(e 3/2 ). (51) 

Noting as well that the averages of A and dA/dtp vanish by construction of 
the transformation, one can average the new system, drop terms of order e 3 / 2 
and higher, and arrive at the second order partially averaged system: 

D = -e 1/2 F ll - eDF 12 , (52) 

G = -eF 22 , (53) 

V=-£ 1/2 -^D + e(-^D 2 + F 32 ), (54) 

g = eF i2 . (55) 

4.2 COMPUTATION OF AVERAGED EQUATIONS 

Computation of the averages of Fij according to equation (43) results in the 
second order partially averaged system: 

D = _ e i/2[^ L 4 (Am + Bm) sin {2g + m(p) + A^ (8 + ™ e 2 + 37 e4)] 



- Jz [L4(j4m + Bm)]lL * sin {2g + mcp) + 3lW {U6 + 37e2)} ' (56) 

G = -e[^Li(A m + B m ) sin (2g + m<p) + -^(8 + 7e 2 )], (57) 



>- -^h D+ * { h D2 -l Ll+ \ L * G2 



3 d 

- g [L\A m + B m )]\ L , cos (2g + m<p)}, (58) 

g = e[- A GLl - + B m ) cos (2.g + rmp)], (59) 

where e = (1 — G 2 /L 2 ) 1 / 2 . Recall that A„ and B v are functions of the eccentric- 
ity e. At resonance the value of L is L* and the type of resonance is governed 
by the resonance relation mL~ 3 = ntt' . Through the integration to determine 
the averages of F^ in equation (43), it is found that the only resonances that 
contribute in the averaging process are the (m : 1) resonances. 
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One can make a simplification by noticing that the arguments of the trigono- 
metric functions all combine as 2g + mip. One can then define a variable $ such 
that <& = 2g + mip. If one multiplies the differential equations for g and ip by 
the appropriate factor and adds them, one has a reduced system of differential 
equations that consists of three equations. One of the variables is now an angle 
$. The system reduces to 

D = -^[yl^il™ + B m ) sin $ + AT a ] 

~ eD{ ^M [Ij4{Am + Bm)]lL - sin<i> + Ar>} ' (60) 

G = -e[^Li(A m + B m ) sin® + AT G ], (61) 

3 T7i d 3 d 

- —qZ^ 4 ^ + b ™)]\l, cos$ - -Lt — (A m + S m )cos$}, (62) 

where the gravitational damping terms are: 

r^^(8+fe 2 + | e 4 ), (63) 

TfJ= 3W (146 + 37e2) ' (64) 
r G = Z ^ I (8 + 7e 2 ). (65) 

The result of the coordinate transformations and averaging is a system of 
three differential equations which describe the averaged motion to second order. 
The resultant set of ODEs represents the evolution of these 'slow' variables. 



5 THE AVERAGED SYSTEM 

The averaged system derived in the previous section is in the form of an ex- 
pansion in the powers of the small parameter e 1 / 2 to second order. One can 
investigate the two different orders of the partial averaging that are calculated 
here, and arrive at results that can be compared to those found numerically. 
The different orders bring to light different features of the solution. 
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5.1 FIRST-ORDER SYSTEM 



To extract the first-order partially averaged system one needs to only keep 
terms of order e 1 / 2 , the small parameter of the perturbation expansion around 
the resonance manifold. One obtains the following system: 

t> = -e^ 2 [^Li(A m + B m ) sin $ + AT Q ], (66) 

G = 0, (67) 
$ = -eV*3m A (6g) 

where A m and B m follow from equations (16) and (17) and T a is defined in 
equation (63). 

One immediate consequence of this system, as can be seen in (67), is that the 
orbital angular momentum G is constant at resonance; hence, the eccentricity 
of the orbit is constant as well. As a result, one can rewrite the system in an 
equivalent way in terms of two canonically conjugate variables D and The 
corresponding Hamiltonian is: 

H = ^D 2 + eAcos$ - er$, (69) 
where A and t are constants 

q m 2 

\ = l—(A m + B m ), (70) 

r = f?AIV (71) 
The canonical equations of the equivalent system are: 

D = eX sin $ + er, (72) 

$ = D. (73) 

Two pendulum-like equations of motion result, one in D and one in The one 
in D traces the orbit of the deviation from the resonant value L* as defined in 
(33) and is expressed as: 

D + eu 2 D D = 0, (74) 

where 

t^ = -Acos$. (75) 

The differential equation in $ describes a pendulum with constant torque. Of 
course, to have oscillatory behavior in both cases, the condition lo 2 d > must be 
met. The first-order approximation of the resonance behavior in D consists of an 
oscillation with slowly varying frequency and constant amplitude. A comparison 
between the analytical formula, particularly that of the frequency (75), and the 
corresponding numerical result can be made to test the quantitative agreement 
between the two methods. 
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5.2 SECOND-ORDER SYSTEM 



When one considers the second-order averaged equations, one can derive a for- 
mula for the damped or antidamped oscillator. It goes as follows: 

D + e'yD + eu 2 D D = 0. (76) 

There is also a corresponding equation for $. The frequency is the same as 
in the previous case, but the coefficient to the damping term can be found from 

7 = l m JZ ^ 4{Am + Bm) ^- Sin<& + A 3gW 146 + 3?e2) ' (?7) 

This describes an oscillation whose amplitude changes with time due primarily 
to damping or antidamping. Also, G is no longer constant in the second order 
as shown in equation (61). While the orbit is in resonance, its orbital angular 
momentum and eccentricity change slowly over time. 

For the purpose of finding an approximate solution to the differential equa- 
tion (76) one can assume that 7 and lojj are constant. That solution is: 

D(t) = Aer^ 1 cos (e 1/2 uj D t + constant). (78) 

The attention will be focused here on D because the numerical data to be 
compared with the averaged results involve the deviations from the resonance 
value, which is what D measures. From this equation, one can see that there is 
a sinusoidal oscillation multiplied by a time-dependent exponential amplitude. 
The character of this function depends on the magnitude and sign of 7 given in 
equation (77). The amplitude increases if 7 < and decreases if 7 > 0. The 
magnitude of 7 affects the rate at which the amplitude changes. More generally, 
7 itself is a function of time; this may result in local variations of the sign of 7. 

6 NUMERICAL MEASUREMENTS 

The measured frequency Lo m is simply an average frequency determined over 
several periods of the numerical L versus t graph. One must keep in mind that 
the averaging is valid over a time interval of eT 1 / 2 . To measure the ratio between 
the value of the envelope of D at two points, that is R m , one simply measures 
the respective values at different times and divides them. 

How does one compare the numerical results to the analytical? First of 
all, it is necessary to recall the derivation of the averaged equations. They are 
derived under the assumption that the orbit is in resonance. Since L = L* 
determines the resonance manifold, the points (L, G, £, g) such that L = L* are 
on the resonance manifold. In the L versus t graph the orbit passes through the 
resonance value as it oscillates. It is these points of the numerical solution 
that are on the resonance manifold that apply to the analytical formulae for ojd 
and 7. 

Recall that the first order averaged equations give a formula for the frequency 
of the solution. It is lo 2 d = —ejK m (e) cos$. This formula depends on e and $ 
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at the resonance manifold. To produce a value of ujd one needs the values of e 
and $ on the resonance manifold. In essence, what is being investigated is the 
agreement between the measured frequency of a resonance recorded in a graph 
and a prediction made by a derived formula. The derived formula is based on 
the orbit in resonance. For the purposes of testing this agreement one can take 
the values of (L, G, £, g) that occur at resonance in the numerical case and apply 
them to the derived formula. 

The solution to the differential equation derived from the second-order case 
gives rise to a formula that describes the behavior of the envelope of the reso- 
nance orbit. To find the analytical formula for this ratio let V(t) be a function 
that describes the envelope of D(t). That is, V(t) = ,4cxp(— \eyt). One can 
compute the ratio between V(ti) and V(t 2 ), where At = t 2 — ii, by dividing the 
two expressions for V(ii) and Vfo)- One comes up with Rd = exp(— ^e^At). 
The determination of 7 calls for the numerical resonance manifold values of e, 
G, and 

The examples that were computed for this study in figures 1-3 yielded the 
following results where Lo m is the angular frequency and R m is the ratio of the 
amplitudes of the envelope at two different times for the (m : 1) resonance. 
Also, the following definitions are used: 

Aw _ \uj d - u m \ 

LO UJD 

AR _ \Rd — R m \ 

R Rd 

The results are: 

Resonance Auj/uj AR/R 

(1 : 1) 0.075 0.088 

(2 : 1) 0.0002 0.0186 

(3 : 1) 0.57 

A condition for the averaging theorem to be valid is that e be sufficiently 
small. Therefore, when computing the numerical results in search of resonances, 
one can pose the question of whether the numerical choice of e is small enough. 
How small is sufficiently small? By lowering the value of e, which in turn leads 
to a more cumbersome and time-consuming numerical calculation, one can most 
likely improve the numerical results. As can be seen in the summary of results, 
the (1:1) and (2:1) examples show good agreement between the numerical and 
analytical results. The (3:1) case, however, deviates quite far from agreement. 
It is known from numerical evidence that the higher-order resonances yield 
more elaborate structure and are more likely to exhibit chaos than low-lying 
resonances (Chicone et al.,1997). In this case, the numerical results could be 
precarious in their agreement with analytical predictions. Perhaps the choice 
for e was not sufficiently small to meet the conditions of the averaging theorem. 
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Figure 4: This is an example of a phase cylinder representing a pendulum 
equation with constant external torque. Specifically, the system is D — asin$ + 
b and $ = D, where a = 1 and b — 0.9. In this example and 2tt are identified 
on the horizontal $ axis. The vertical axis is D = $. The elliptical orbits 
represent periodic solutions around a rest point. There is a seperatrix beyond 
which the solutions are not periodic. 

7 ANALYSIS OF RESONANCE CAPTURE 

Librational motion on the phase plane of the first-order model is indicative of 
an orbit captured into resonance (Chicone et al., 1997). The pendulum- with- 
torque system shown in (72) and (73) must have an elliptical fixed point to 
have these librations. Furthermore, a hyperbolic fixed point is expected with 
a homoclinic orbit to give rise to an area in the phase space where orbits pass 
through resonance. A fixed point on the phase plane is given by (D, $) = (0, $ ), 
where $ satisfies: 

sin$ = -T/A. (79) 

In general, three different scenarios can arise as the result of the pendulum 
system. For there to be a solution for $o the inequality |r/A| < 1 must hold. 
For the cases where |r/A| < 1 or t/A = 1, there are fixed points. Consider the 
case where the inequality |r/A| < 1 holds. The phase cylinder for this shows 
that there exist orbits that are oscillatory and stay in the resonance capture 
region, and also orbits that pass through the resonance (see Figure 4). This 
shows that whether the orbit will be captured into resonance depends on the 
initial conditions. The case such that |r/A| > 1 can arise where there is no 
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solution to the equation and therefore no fixed point. 

The first-order system of equations that gives rise to these pendulum- like 
equations displays in an approximate way the physics of the original astrophys- 
ical model. From the first-order system one can see from equation (71) that 
the torque r is linearly related to the gravitational radiation damping by the 
coefficient A. If there is sufficiently strong damping, and hence a large A, there 
will be no resonance capture because r/A can in all cases exceed unity. All 
solutions will pass through the resonances and the orbit will follow a path such 
that the semimajor axis will shrink to zero. This would be the expected result 
if gravitational radiation damping were the sole perturbation. Recall also that 
the binary system will be captured only into the (m : 1) resonances according 
to the averaged equations. 

The second-order averaged system includes antidamping effects as well as a 
slowly drifting A and r (Chicone et al. 2000). It offers a description of how 
the actual system can fall out of resonance by leaving the capture region by its 
anti-damped motion. The full non-linear system exhibits effects like resonance 
capture and exit from resonance, dynamics whose full description elude current 
analytical techniques. However, the second-order system explains with some 
success the system's behavior in resonance, particularly the frequency and the 
rate of antidamping. 

There are interesting features in the (3 : 1) resonance, particularly its alter- 
nating sets of grouped oscillations in an 'excited' phase and 'relaxed' phase. The 
slowly changing variables of L, G, and g, and the fast variable £ as seen in the 
Delaunay form of the system (6)- (9) make it a candidate for the phenomenon 
of bursting. One may speculate from the characteristic features in figure 3 that 
bursting is present in this (3 : 1) resonance (Izhikevich 2001). 

The phenomenon of resonance capture describes the system's state where the 
dissipative influence of gravitational radiation damping is offset on the average 
by the tidal perturbation of the orbiting third mass; that is, the theory of 
capture into resonance describes how the orbit's loss of energy via gravitational 
radiation reaction is countered by a deposit of energy from the orbit of the third 
mass. 

8 CONCLUSION 

The averaging method produces a set of ordinary differential equations that 
describe the behavior of the system while in resonance. The quantity D tracks 
the oscillatory motion of L about L*. The level of approximation given by 
the first-order averaged equations highlights the oscillatory characteristics of 
the actual solution whereas the second-order averaged equations reveals the 
slow changes in amplitude. The numerical solution shows more structure and 
detail that is inherent in the actual solution. There are prominent features 
of the actual solution that one can compare qualitatively and quantitatively 
to that of the averaged system. Though the first-order solution only predicts 
oscillations of constant amplitude and slowly changing frequency and the second- 
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order solution predicts an exponential trend of damping or antidamping, they 
correspond on the appropriate timescale to the most prominent features of the 
numerical solution. This makes sense in the light of what the averaging method 
intends to do — capture the 'slow' motion. 

Resonance capture and chaos are two signature results that arise from a 
system in resonance. Perhaps the most noted example of this phenomenon is 
the damped driven oscillator. The astrophysical system that is modelled in this 
study by a damped and driven Kepler system displays similar characteristic 
effects — namely resonance capture. Under the non-Hamiltonian perturbation, 
the nonlinear Hamiltonian system can exhibit both resonance capture and chaos 
(Chicone et al. 2000). This can be manifest in the attempt to compare the 
numerical case with the analytical. An orbit near a higher-order resonance such 
as (3 : 1) will more likely be chaotic and therefore prone to bear 'erroneus' results 
when compared with the partially averaged equations. The dense structure 
shown in figure 3 that corresponds to (3 : 1) suggests such chaotic motion. 
However, good results followed from the analysis of the (1 : 1) and (2 : 1) 
resonances. 

The averaged system to first order shows that at each (m : 1) resonance 
one sees motion like a pendulum with constant torque. To the second order, 
antidamping and nonlincarity may cause disruption to the resonance. The full 
nonlinear system includes more effects that are not seen in the averaged sys- 
tems which may contribute to the disruption of the resonance. The analytical 
formulas for the first and second-order averaged systems outline structure that 
is part of the actual solution. It is also of interest to determine how well the 
results from the analytical part agree with the numerical. The complexity of the 
actual nonlinear equations and their solution makes analytical formulae helpful 
and illuminating. Both analytical and numerical techniques given here present a 
consistent description of the nonlinear effects associated with resonance capture. 
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